library(sp)
library(tidyverse)
library(readxl)
library(rgeos)
library('raster')
library('geosphere')
library(maptools)
library(maps)
library(haven)
options(repos=c(CRAN="http://mirror1.csvd.census.gov/CRAN"))
#install.packages('rgdal')
require ('rgdal')

# this data is each turbine, not just each plant: 
setwd("//projects//")
uswindplants <- read_dta("users//########//PlantLevelWindDatabaseCSM.dta")

# number of meters in each d-mile radius (don't really need this if run above)
maxdist_5 <- 8047 
maxdist_10 <- 16093.5
maxdist_20 <- 32187    
maxdist_40 <- 64374
maxdist_60 <- 96561
maxdist_80 <- 128748
maxdist_100 <- 160935

turb<-data.frame(cbind(uswindplants$xlong, uswindplants$ylat)) 
names(turb)<-c("Lon_T","Lat_T") 

# convert wind turbines to spatial object
t.spoint.p <- SpatialPointsDataFrame(turb[, c("Lon_T", "Lat_T")], 
                                     data.frame(ID=seq(1:nrow(turb))), 
                                     proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
# transform to planar projection: 
t.spoint.mp <- spTransform(t.spoint.p,CRS( "+init=epsg:2163" )) 

for (k in c(20,40,60,80,100)) {
  #tbuf.k.p <- gBuffer(t.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE)
  assign(paste0("tbuf.",k,".p"),gBuffer(t.spoint.mp, width=get(paste0("maxdist_",k)), byid=TRUE))
  assign(paste0("t.",k,".p"), gUnaryUnion(get(paste0("tbuf.",k,".p")), id = NULL))
  assign(paste0("t.",k,".p.reproj"), spTransform(get(paste0("t.",k,".p")), CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
  assign(paste0("t.",k,".p.reproj"), as(get(paste0("t.",k,".p.reproj")), "SpatialPolygonsDataFrame"))
  # Subset to household locations within each of these buffer ranges (not sure if can do this with census data)
  # assign(paste0("hh_in_",k),h.spoint.p[get(paste0("t.",k,".p.reproj")),])
}

##Shapefiles for larger k-mile buffers##
writeOGR(t.20.p.reproj, dsn= "/projects/programs", layer= '20mileUSp', driver= 'ESRI Shapefile')
writeOGR(t.40.p.reproj, dsn= "/projects/programs", layer= '40mileUSp', driver= 'ESRI Shapefile')
writeOGR(t.60.p.reproj, dsn= "/projects/programs", layer= '60mileUSp', driver= 'ESRI Shapefile')
writeOGR(t.80.p.reproj, dsn= "/projects/programs", layer= '80mileUSp', driver= 'ESRI Shapefile')
writeOGR(t.100.p.reproj, dsn= "/projects/programs", layer= '100mileUSp', driver= 'ESRI Shapefile')

#map('state')
#map(SpatialPolygons2map(t.100.p.reproj),add=TRUE,col=2)
